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Using a Monte Carlo (MC) method, we study an effective model for the Fe(II)Fe(III) bimetallic 
oxalates. Within a hybrid quantum-classical MC algorithm, the Heisenberg S = 2 and S' = 5/2 
spins on the Fe(II) and Fe(III) sites are updated using a quantum MC loop while the Ising-likc 
orbital angular momenta on the Fe(II) sites are updated using a single-spin classical MC flip. The 
effective field acting on the orbital angular momenta depends on the quantum state of the system. 
We find that the mean-field phase diagram for the model is surprisingly robust with respect to 
fluctuations. In particular, the region displaying two compensation points shifts and shrinks but 
remains finite. 

PACS numbers: 75.50.Xx, 71.70.Ej, 75.10.Dg, 75.40.Mg 
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I. INTRODUCTION 

Bimetallic oxalates are layered, molecule-based mag- 
nets with the chemical formula A[M(II)M'(III)(ox) 3 ] [if. 
Every layer contains two different transition metal atoms, 
M(II) and M'(III), in the alternating honeycomb struc- 
ture depicted in Fig[TJ Each bond represents an oxalate 
molecule ox = C2O4, which generates a crystal-field po- 
tential at both ionic sites. For different transition metals, 
bimetallic oxalates can be ferromagnetic, antiferromag- 
netic, or ferrimagnetic with moments always pointing out 
of the plane @. Since the type of magnetic order does 
not depend on the cation A that couples the magnetic 
layers, the magnetic properties of the bimetallic oxalates 
are primarily controlled by a single bimetallic layer. 

For the Fe(II)Fe(III) bimetallic oxalates, however, the 
presence of magnetic compensation below the ferrimag- 
netic transition temperature T c does depend on the 
choice of A 0]. For several cations, the magnetization 
in a small field is positive just below T c as 45 K but then 
become negative below r comp ss 32 K. This effect was 
explained by Fishman and Reboredo [J] , who used mean- 
field (MF) theory to solve an effective Hamiltonian that 
includes spin-orbit coupling on the Fe(II) sites. Magnetic 
compensation is produced when the orbital angular mo- 
mentum L cf of the low-lying crystal-field doublet on the 
Fe(II) sites exceeds a threshold value. By altering the 
crystal-field potential, the cation A can shift L cf above 
or below this threshold value. In this paper, we use a 
Monte-Carlo (MC) technique to study the same effective 
Hamiltonian and demonstrate that the MF results are 
surprisingly immune to the effect of fluctuations. 

By Hund's first rule, the spin on the Fe(II) (3d 6 )and 
Fe(III) (3d 5 ) sites are S = 2 and S' = 5/2, respectively. 
Since the Fe(III) multiplet is half-full, its orbital angu- 
lar momentum L' vanishes according to Hund's second 
rule. The L — 2 orbital angular momentum on the Fe(II) 
sites is split by the GVsymmetric crystal-field potential 
produced by the 6 oxygen atoms surrounding each ion. 



This splitting creates two doublets and one singlet [4(. 
The orbital angular momentum is unquenched when one 
of the doublets lies lowest in energy. In that case, the 
out-of-plane or z component of the orbital angular mo- 
mentum on the Fe(II) sites takes values ±L cf , where L ci 
ranges from to 2 and depends on the crystal-field po- 
tential. The spin-orbit coupling on the Fe(II) sites is 
given by AL; • Si, where L, = ±L cf z and A w —147 K 
[5j is the spin-orbit coupling constant (negative because 
the 3d 6 shell is more than half-filled). The total angular 
momentum J; = Si + Lj on the Fe(II) sites is not a good 
quantum number and Hund's third rule is not obeyed 
because the crystal-field potential is large compared with 
the spin-orbit coupling. 

Hence, the effective Hamiltonian of the Fe(II)Fe(III) 
bimetallic oxalates can be written 



•S' • A^.V'. 



(1) 



where the (ij) summation in the exchange term is per- 
formed over all nearest neighbors on the honeycomb lat- 
tice and the i summation in the spin-orbit term is per- 
formed over the Fe(II) sites only. The orbital angular 
momentum h\ — ±L ci can be treated as a classical vari- 
able so long as the relevant energy scales are smaller than 
the splitting A between the lowest-energy doublet and 
the nearest excited states of the crystal-field potential. 
Since T c is less than 45 K while A is larger than room 
temperature, this should be a good assumption for the 
bimetallic oxalates. For the special case where the singlet 
lies lowest in energy, we would take L cl = 0. 

Several other models can also explain the existence of 
magnetic compensation in the Fe(II)Fe(III) bimetallic ox- 
alates. Nakamura [6| described the Fe(II) and Fe(III) 
spins by S = 2 and S' = 5/2 Ising variables with single- 
ion anisotropy —D(Sf) 2 on the S — 2 sites. Although 
magnetic compensation is found in a MF or effective- 
field theory treatment of this model, it is absent when 
the model is solved using the MC method. With flue- 
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FIG. 1: The honeycomb lattice showing the alternating Fe(II) 
and Fe(III) sites. 



tuations correctly included, magnetic compensation is 
recovered only after interlayer or longer-ranged interac- 
tions are considered. A similar conclusion was reached 
by Carling and Day Q, who used MC simulations to 
show that next-neighbor interactions between Ising spins 
on the same sublattice are required to obtain magnetic 
compensation within a single bimetallic layer. Li et al. 
[H treated the spins as Heisenberg operators with inter- 
layer interactions and single-ion anisotropy on both sub- 
lattices. Using a Green's function decoupling scheme, 
they found that magnetic compensation occurs even in 
the absence of interlayer interactions when anisotropy is 
experienced by a single magnetic sublattice. However, it 
is uncertain whether magnetic compensation would sur- 
vive in the presence of correlated thermal fluctuations. 

Those other treatments fail to explain several impor- 
tant features of the Fe(II)Fe(III) bimetallic oxalates that 
can be explained by the model of Eq.JT]). First, Eq.JT]) 
provides a natural explanation for the origin of the mag- 
netic anisotropy on the Fe(II) sites due to the splitting 
of the L = 2 multiplet by the crystal-held potential. 
Second, unlike the model studied by Nakamura @, the 
model studied in this paper does not rely on interlayer in- 
teractions to produce magnetic compensation. So it can 
explain the appearance of magnetic compensation in the 
A = N(n-C„H2„+i)4 family: as n increases from 3 to 5, 
the interlayer separation grows from 8.2 A to 10.2 A Q 
but magnetic compensation appears only for n = 4 and 
5, probably because defects in the bimetallic layer are 
created by the shortest cation. Third, Eq.fJJ can explain 
the persistence of negative magnetization below T comp in 
small holds: flipping the orbital angular momentum L\ 
requires an energy of 2\X\L cl S « 175 K ^ T comp . By con- 
trast, the single-ion anisotropy D ~ A 2 /A experienced 
by the Fe(II) moments in most materials is less than 10 



K 0. So the energy barrier AD for flipping the Fe(II) 
spin is smaller than about 40 K, which is comparable 
to T comp ~ 32 K. Finally, Eq.fJJ explains the recently- 
observed jump (lo| in the magnetization between T comp 
and T c , which is believed to arise from an inverse Jahn- 
Teller transition [llj. Indeed, it is difficult to explain this 
jump based on the other models introduced above. 

Magnetic compensation occurs in a MF treatment [i[ 
of Eq.([T]) for two reasons: the anisotropy at the Fe(II) 
sites produced by spin-orbit coupling and the Fe(II) or- 
bital contribution to the total magnetic moment. The 
latter is absent in an Ising or Heisenberg ferromagnet 
with single-ion anisotropy. While MF theory was used 
to approximate the exchange coupling JSi ■ Sj between 
neighboring Fe(II) and Fe(III) moments, the spin-orbit 
coupling XLi ■ Sj at the Fe(II) sites was treated exactly 
within the crystal-field doublet. As demonstrated by 
Nakamura Q, however, models that exhibit magnetic 
compensation within MF or effective-field theories may 
no longer do so once fluctuations are correctly included. 
So a MC study including both quantum and thermal fluc- 
tuations is needed to confirm that Eq.JTJ supports mag- 
netic compensation. 

This paper is divided into four sections. In the next 
section, we describe the hybrid MC method used to study 
Eq.([T]). The results of our MC study and a comparison 
with earlier MF results are presented in Section III. A 
brief conclusion is provided in Section IV. 

II. HYBRID MC METHOD 

We apply the stochastic series expansion (SSE) quan- 
tum MC method [l^] to the model Hamiltonian H. The 
SSE method employs a Taylor expansion of the partition 
function Z: 

^EE^H^H, (2) 

a n—0 

where \a) are the basis states used to evaluate the ma- 
trix elements of H and (3 — 1/T is the inverse temper- 
ature. Quantum MC methods have been previously ap- 
plied to high-spin models and here we follow the method 
described in Ref.[13j. As explained below, that method 
has been modified to treat the Ising-like orbital angular 
momentum. For a more detailed account of the general 
SSE method, we refer to Ref-Q. 

To formulate the updating procedure, we write the 
Hamiltonian as a sum over all nearest neighbors in the 
system 

H=-jJ2H ij - (3) 

The operator Hij can be decomposed into its diagonal 
and off-diagonal parts: 

Hij = Hn^ij — Ho,ij , (4) 
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where the subscript D denotes a diagonal operator and 
O an off-diagonal operator. For the present model, these 
two operators take the form 



Hdai — C — SfSj 



A 

3J" 



S^ L~ 



and 



Ho,ij — 



5^ 



s, 



(•5) 



(6) 



The constant C is included in order to ensure a positive 
weight in the expansion. 

Introducing a cutoff K in the Taylor expansion (which, 
when done properly, does not cause any systematic errors 
[Til ]) and including additional unit operators /, the ex- 
pansion can be rewritten as 



(3 n (K -n)\ 
~K\ 



(a\SK\a), 



(7) 



where Sk is the operator string 

K 

Sk = Y\_ Hp> 

p=i 



(8) 



with Hp G {Hn,ij, Ho,ij, I}- Now n is the number of 
bond-operators Hjj^j or Ho,ij in the operator string Sk 
[ijj]. The MC procedure must sample the space of all 
states \a) and all operator sequences Sk with the relative 
weight 



W(a, Sk) = (a\S K \a) 

Denoting a propagated state by 



(9) 



while the probability for removing a diagonal operator is 



Pr, 



K - n + 1 



N p f3{a{p)\H D>ij \a(p)y 



(12) 



The total number of nearest-neighbor pairs on the lat- 
tice is denoted N p . In a diagonal update, one exchange 
attempt is made for each diagonal and unit operator. 

The second type of update is a global operator-loop 
update, which leaves unit operators unaffected. This up- 
date forms and flips a closed loop of spins in the vertex 
list. In this process, both the affected vertices and states 
are changed. For a detailed description of the operator- 
loop move, we refer to Ref. [lj|. In the absence of spin- 
orbit coupling, the operator-loop update together with 
the above diagonal update ensure that the complete SSE 
configuration space is sampled. But when spin-orbit cou- 
pling is included, the Ising-like orbital angular momen- 
tum Lf must also be updated. Due to the classical na- 
ture of Lf , no terms in the Hamiltonian are able to flip 
Lf and the orbital angular momenta must be updated in 
a separate move, which we describe next. 

Additional flips of the spin variables S* and S'j z are 
allowed in case no string operator acts on sites i or j. 
The weight of the configuration then remains unchanged 
after a spin flip. The same is true for the orbital angular 
momenta Lf : if no string operator acts on a given vari- 
able Lf , it can be flipped with no associated change in 
weight. 

This method must be modified at or below T c , where 
many string operators act on a given Lf. Flipping Lf 
then causes the corresponding vertices to change with the 
associated weight change given by Eq. , which must be 
must be taken into account when updating the orbital 
angular momenta. Denoting the weight of a vertex by 



W p = (a(p)\H p \a(p-l)), 



(13) 



Hp)) = n^l a >< 



(10) 



1=1 



the matrix element in Eq. ([7]) can be written as a product 
of elements with the form (a(p)\Hij\a(p — 1)), which is 
equivalent to 



(S?(p)Lt (p)S'f (p)|ff tf \S?(p - l)L?(p-l 



1)>- 



We shall refer to these matrix elements as "vertices." The 
matrix element in Eq.© can be viewed as a list of such 
vertices. 

In the operator-loop algorithm, two basic updates en- 
sure that the complete SSE space is sampled. The di- 
agonal update attempts to exchange diagonal operators 
Ho,ij with unit operators /. The probability for inserting 
a diagonal operator (exchanging it for a unit operator) 
at position p in the operator sequence is 



N p [3{a(p)\H D ^\a{p)) 
K -n 



(11) 



then the acceptance probability of flipping an orbital mo- 
mentum is 



P = 



(i) 



(14) 



p(i) 



where W indicates the weight after the flip and W the 
weight before the flip. The product runs over all vertices 
containing the orbital momentum h\ . 

The method described here works well over the range 
of parameters studied in this paper. It acts like a hy- 
brid quantum-classical MC method where the Heisenberg 
spins are updated using a quantum MC loop algorithm 
but the Ising-like orbital angular momenta are updated 
separately using a classical single spin-flip MC method. 
The effective field experienced by each orbital angular 
momentum depends on the quantum state {\o),Sk} of 
the system. Since the orbital angular momenta saturate 
before the Heisenberg spins, it may be preferable to turn 
off the classical update at low temperatures in order to 
prevent freezing in a state different from the ground state. 
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It is straightforward to generalize the present method to 
include orbital angular momenta on both magnetic sub- 
lattices [HI , which would allow the study of many other 
families of bimetallic oxalates. Like the standard stochas- 
tic series method, this method can also easily be gener- 
alized to arbitrary spin size. 



III. TRANSITION TEMPERATURE AND 
PHASE DIAGRAM 

We begin our discussion by considering the lattice and 
boundary conditions. The honeycomb lattice consists of 
two sites per unit cell. Using periodic boundary condi- 
tions on a two-dimensional lattice with N x N unit cells, 
the number of sites iV s = 2N x N varies between 32 
(N = 4) and 32768 (N = 128). To study much larger 
system sizes would most likely require improvements in 
the algorithm. 

The thermodynamic expectation value of the sublat- 
tice magnetization on the Fe(II) sites is given by 



(\s z \) = 



(15) 



where the summation runs over all N s /2 sites. In a sim- 
ilar manner, we define the sublattice magnetization on 
the Fe(III) sites, (|5" z |), and the orbital angular momen- 
tum at the Fe(II) sites, (|£ 2 |). The absolute value is 
measured because the strict statistical average vanishes 
for finite systems due to time-reversal symmetry in zero 
magnetic field. The total magnetization is given by 



M = 2(\S'*\)- 2<n> -<|L«|>, 



(16) 



where we take 3 = 2. Recall that (S z ) and (L z ) have the 
same sign because A is negative. 

Sublattice magnetization curves for A = —13.3 J and 
L ci = 0.3 are plotted versus temperature in Figj2] As 
shown later, these parameter values should be close to 
those expected for bimetallic oxalates that exhibit mag- 
netic compensation. Many features characteristic of the 
general magnetization curves appear in this figure. At 
low temperatures, the results for different system sizes 
have converged to the thermodynamic limit. But due 
to the diverging correlation length, increasing finite-size 
effects prevent an accurate determination of the critical 
temperature. For the remainder of this paper, we only 
present those portions of the magnetization curves that 
have converged in system size. 

At low temperatures, we can clearly see the effects 
of quantum fluctuations. Since the spins are treated as 
Heisenberg operators, the classical ground state is not an 
eigenstate of the Hamiltonian. The ground state magne- 
tizations for the S = 2 and S' = 5/2 sublattices are 
1.91 and 2.41, respectively, or about 4% below the clas- 
sical values. These deviations are small because of the 
relatively high spin values. The figure also illustrates 
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FIG. 2: Sublattice magnetizations (\S'\), (\S\) and (\L Z \) as a 
function of temperature for A = —13. 3 J and L ct — 0.3. The 
solid and dashed lines denote results for 4x4 and 64 x 64 unit 
cells, respectively. The critical temperature is T c = 4.28J 



the physical mechanism responsible for magnetic com- 
pensation. Due to the effective anisotropy induced by 
the spin-orbit coupling to the Ising-like orbital angular 
momentum, the S = 2 sublattice magnetizes faster than 
the S' = 5/2 sublattice. At T = 2 K, the S sublattice 
has almost reached its saturation magnetization but the 
S' sublattice has not. 

In order to accurately determine the critical tempera- 
ture, we have measured the Binder ratio [l6| of the mag- 
netization, 



Q = 



(17) 



which is the ratio of two moments of the order parame- 
ter. The Binder ratio should become size-independent at 
a second-order phase transition. After plotting the ratio 
for various system sizes, the point of intersection gives 
the critical temperature. In Fig[3J we plot the Binder 
ration for AL cf = — 4J, which is the value used in FigfS] 
Finite size effects are clearly visible: the ratios for sys- 
tem sizes 4x4 and 8x8 intersect at T = 4.34J while 
the ratios for larger system sizes intersect at T = 4.28J. 
The convergence in system size improves with increasing 
anisotropy. 

The critical temperature depends only on the prod- 
uct — XL ci and we compare the MF and MC results for 
T c in FigJU According to the Mermin- Wagner theorem 
[13], the critical temperature of a two-dimensional sys- 
tem must approach zero as the anisotropy vanishes. One 
of the main weaknesses of the MF solution plotted in 
FigfJ] is that it overestimates T c and predicts a finite 
critical temperature even for an isotropic system with 
— XL = 0. By contrast, the MC results in Figf?]indicate 
that T c rises very rapidly with small —\L ci . However, 
even a small anisotropy induces a critical temperature of 
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FIG. 3: Binder ratio as a function of temperature for \L ci = 
— 4J. Results for different system sizes intersect at the critical 
temperature T c = 4.28J The number of sites varies from 32 
(N = 4) to 8192 (N = 64). 
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FIG. 5: Total and sublattice magnetization calculated in MF 
theory (dashed line) and in the MC method (solid line) for 
A = -12 J and L ci = 0.25. We display only MC data that 
have converged in system size, therefore not quite reaching 
the critical temperature. 
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FIG. 4: Critical temperature as a function of spin-orbit cou- 
pling — \L ci . Results from MF theory and quantum MC 
calculations along with the asymptotic behavior for small 
anisotropy from RG theory. 



order J. In FigHJ we model the rapid initial rise in the 
critical temperature with the functional form 



HWIMY 



(18) 



This dependency was first reported by a renormalization- 
group (RG) study [18|, where the parameters a and b 
were related to the critical temperature of the three- 
dimensional isotropic model and the the anisotropy pa- 
rameter, respectively. Here, we simply treat a and b as 
fitting parameters to demonstrate that our MC results 
are reasonable. Notice that MF theory overestimates T c 
by about 60% in almost the whole parameter range of 
interest, an effect that must be taken into account when 
applying MF theory. 



Scaling the temperature by T c , we compare the MC 
and MF results for the total and sublattice magnetiza- 
tions in FigE for A = —12 J and L ci = 0.25. Below 
0.5T C , the classical variable L z is essentially in its fully- 
polarized classical ground state within both the MF and 
MC solutions. For the Heisenberg variables S and S', 
quantum fluctuations again cause deviations from the 
classical ground state. In the critical region, the sub- 
lattice magnetizations are described by the functional 
form (T - T C )P. Whereas the MF exponent is (3 = 1/2, 
the Ising-like exponent (3 — 1/8 describes the MC data. 
Consequently, MF theory underestimates both sublattice 
magnetizations as T c is approached. The agreement be- 
tween the MF and MC results for the total magnetization 
M = 2(\S' Z \)-2(\S Z \)-(\L Z \) is generally better than for 
the individual sublattice magnetizations. Since the total 
spin 5 t z ot = Yli + J2j S'f commutes with H, quan- 
tum fluctuations do not effect the total magnetization at 
T = 19] and M(0) = 2(5" - S) - L ci . At higher tem- 
peratures, the MF solution underestimates M(T), which 
has the same critical exponent (3 = 1/8 as the sublattice 
magnetizations. 

Due to both the orbital contribution and the effects of 
anisotropy, the total magnetization depends sensitively 
on L cl , as seen in FigJB] for A = — 8J. When L cl is 
small, M is dominated by the S' = 2 sublattice (defined 
to be positive) and it remains positive for all temper- 
atures. When L cf is large, M is dominated by the S 
sublattice and it remains negative for all temperatures. 
For intermediate values of L cf , there is competition be- 
tween the two sublattices. As shown, the magnetizations 
for L cf — 0.625 and 0.875 exhibit a compensation point 
Tcomp where M{T) vanishes due to a cancellation on the 
two sublattices. 

By varying the temperature or some other control pa- 
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FIG. 6: Total magnetization for A = — 8J as a function of L c{ 
and temperature. 



rameter such as strain the magnetization can be 
switched in the vicinity of a compensation point. Due 
to the potential applications of this effect, one main goal 
of this study was to determine the values of A and L ci 
where magnetic compensation may occur. Our main re- 
sult is the phase diagram of Fig[7| denoting the number 
n comp of compensation points as a function of L ci and A. 
This phase diagram was previously determined using MF 
theory In FigJTl we display the MC results on top of 
the MF phase diagram. Since the magnetization of the 
classical ground state is given by M = 1 — L ct , L cl = 1 
is an important dividing line: for L cf < l,the S' = 5/2 
sublattice dominates and the ground state magnetiza- 
tion is positive; for L cl > 1, the S = 2 sublattice dom- 
inates and the ground state magnetization is negative. 
When L ct < 1 and —A is large enough, the S sublattice 
may magnetize faster than the S' sublattice, producing a 
compensation point T comp as the temperature is lowered. 
Similarly, when L ci > 1 and —A is sufficiently small, the 
S' sublattice may magnetize faster than the S sublattice, 
again causing the two sublattice magnetizations to cancel 

&f ^corap *^ T c . 

The MF phase diagram contains an interesting region 
with not one but two compensation points. From FigjTl 
it is evident that this area has shrunk and shifted, but 
remains finite despite the effect of fluctuations. It is im- 

= 2 



comp 



possible to say whether the MC region with n 
survives in a narrow neck for all — A/ J greater than about 
1.1. But for —A/ J > 4.5, there no discernible neck of 
L cf where two compensation points can be found. So it 
came as a surprise when a recently-studied Fe(II)Fe(III) 
bimetallic oxalate seemed, at first sight, to exhibit two 
compensation temperatures 10]. However, those mea- 
surements are more naturally explained by an inverse 
Jahn- Teller transition, above which the C3 symmetry of 
the lattice is violated [ill ]. 

For L cf < 1, the MF curve separating the n com p = 2 
and 1 region lies quite close to the MC curve separating 



FIG. 7: Phase diagram indicating the number of compensa- 
tion points n comp . The phase boundaries are given by the 
L cf = 1 vertical line, as well as by the dashed (MF) and 
solid curves (MC). Two compensation points are found in the 
shaded regions. In the inset are shown the MC magnetization 
curves for the three indicated points. 



the n comp = and 1 region when —A/ J is large. These 
two curves appear to cross at about —A/ J « 12. For 
L cf > 1, the difference between the MF and MC regions 
with n comp = 1 is much more pronounced, with the MC 
region about half the size of the MF region. 

In the inset to Fig[TJ we show the magnetization curves 
for three selected points in the phase diagram. These 
points display either one (L cf > 1 or L cl < 1) or two 
compensation points. Notice that the S = 2 sublattice 
dominates at low temperatures for the curve showing one 
compensation point with L ci > 1. 

Using the experimental results, we may now esti- 
mate the parameter values that are relevant for the 
Fe(II)Fe(III) bimetallic oxalates. With A = -147 K [| 
and T c = 45 K, we find that |A|/T C w 3.3. The bimetal- 
lic oxalates should lie quite close to the phase boundary 
in Fig[7] because some exhibit a compensation point but 
others do not 0]. In Fig[2l we used the parameter values 
A = -13.3J and L cf = 0.3, corresponding to T c = 4.24J 
and |A|/T C w 3.13, which is close to the expected ratio. 
Since this point lies near the phase boundary in FigjTl 
it represents a reasonable estimate for the bimetallic ox- 
alates. These estimates yield an exchange constant of 
J w 10.6 K, about twice the MF estimate 0. 



IV. CONCLUSION 

We have studied a model for a molecule-based mag- 
net using a hybrid quantum-classical MC method. Our 
main result is the phase diagram in Figjjj which displays 
the number of compensation points as a function of the 
spin-orbit coupling and orbital angular momentum. The 
region with two compensation points, previously found 
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in MF studies has shrunk but remains finite. 

With a simple scaling of the transition temperature, 
almost all MF results are qualitatively recovered in this 
MC study. The surprising robustness of the MF results to 
the effects of fluctuations arises from two factors. First, 
the spin values S = 2 and 5' = 5/2 in the Fc(II)Fe(III) 
bimetallic oxalates are high enough that quantum fluctu- 
ations are relatively unimportant. Second, the important 
magnetic properties depend on the difference between the 
magnetic moments on the two sublatticcs. That differ- 
ence is more immune to the effect of fluctuations than 
the individual sublattice moments. 

Several of the discrepancies between the earlier MF 
results and experiments are resolved by our MC re- 
sults. Within MF theory, the Curie-Weiss tempera- 
ture O obtained from the high-temperature susceptibility 
X ~ C/(T — 0) is always smaller in magnitude than T r [41 
whereas experimentally it can be twice as large 0, l20j |. 
While 1 6) | should remain close to its MF value even in 
the presence of fluctuations, T c is reduced by about half 
from its MF value. Therefore, the ratio T c /|0| is also 
suppressed by about 50% due to fluctuations. 

The MC results for T c in Figf?] also explain why 
bimetallic oxalates that exhibit magnetic compensation 
tend to have transition temperatures roughly f K higher 
than those that do not [3]. If materials with magnetic 
compensation have — XL cf ~ 4J, then FigfJ] suggests 
that materials without magnetic compensation will have 
— XL ci ~ 2 J with T c reduced by about 25%. By contrast, 
the MF result for T c in FigU] is relatively insensitive to 
— XL and cannot explain this sizeable suppression of the 
transition temperature. 

Since our model Hamiltonian does not contain any in- 
tcrlayer couplings, T c vanishes if the singlet lies below 
the doublets in the L = 2 Fe(II) multiplet. Therefore, 
it seems likely that one of the doublet always lies lowest 
in energy even in compounds that do not exhibit mag- 
netic compensation. Of course, neglecting the interlayer 
coupling was just a convenient approximation within the 
current treatment. While the relative insensitivity of the 
transition temperature to the interlayer separation and to 



the presence of radical spin- 1/2 cations between the lay- 
ers |2(], HH suggests that the interlayer coupling is small, 
it must be present to support long-range magnetic order 
along the out-of-plane direction in zero magnetic field. 

The present method can be easily generalized to con- 
sider spins of any size and spin-orbit coupling on both 
magnetic sublattices. Since quantum fluctuations grow 
with decreasing spin, the hybrid classical-quantum MC 
technique may be particularly useful when considering 
low-spin bimetallic oxalates such as Ni(II)Mn(III) with 
S = 1 and S' = 2. This would enable us to re-evaluate 
the magnetic phase diagrams obtained by Reis et al. [l5| , 
where compensation was found in ranges of L and L' 
but quantum and thermal fluctuations were not consid- 
ered. In order to improve the performance of the method, 
particularly close to T c , one could use the directed- loop 
method L4| to minimize the backtracking process in the 
loop construction. Furthermore, the classical single- spin 
flip update could be augmented by a cluster update [22| • 
We hope that the hybrid quantum-classical MC tech- 
nique developed in this paper also proves of value in 
studies of other systems with both quantum and clas- 
sical degrees of freedom. 
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